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Abstract 

We consider the one- dimensional Maxwell equations with low contrast periodic 
linear refractive index and weak Kerr nonlinearity. In this context, wave packet initial 
conditions with a single carrier frequency excite infinitely many resonances. On large 
but finite time-scales, the coupled evolution of backward and forward waves is governed 
by nonlocal equations of resonant nonlinear geometrical optics. For the special class 
of solutions which are periodic in the fast phase, these equations are equivalent to an 
infinite system of nonlinear coupled mode equations, the so called extended nonlinear 
coupled mode equations, xNLCME. Numerical studies support the existence of long- 
lived spatially localized coherent structures, featuring a slowly varying envelope and a 
train of carrier shocks. Thus, it is natural to study the localized coherent structures of 
xNLCME. 

In this paper we explore, by analytical, asymptotic and numerical methods, the 
existence and properties of spatially localized structures of the xNLCME system, which 
arises for a refractive index profile consisting of periodic array of Dirac delta functions. 

We consider, in particular, the limit of small amplitude solutions with frequencies 
near a band-edge. In this case, stationary xNLCME is well-approximated by an infinite 
system of coupled, stationary, nonlinear Schrodinger equations, the extended nonlinear 
Schrddinger system, xNLS. We embed xNLS in a one-parameter family of equations, 
xNLS e , which interpolates between infinitely many decoupled NLS equations (e = 0) 
and xNLS (e = 1). Using bifurcation methods we show existence of solutions for a range 
of e € (— eo, eo) and, by a numerical continuation method, establish the continuation of 
certain branches all the way to e = 1. Finally, we perform time-dependent simulations 
of truncated xNLCME and find the small-amplitude near-band-edge gap solitons to be 
robust to both numerical errors and the NLS approximation. 
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1 Introduction and Overview 



Nonlinear waves in periodic structures have been a subject of great interest for many years. 
Early interest arose from the possibility of balancing the band dispersion of the periodic 
structure with the nonlinearity to form soliton-like structures; see, for example, [6 10 and 



references cited therein. While such a heterogeneous medium possesses the same soliton- 
producing ingredients of dispersion and nonlinearity as found in the well known Korteweg-de 
Vries (KdV) and nonlinear Schrodinger (NLS) equations which govern nonlinear dispersive 
waves in spatially homogeneous media, the periodic variations of such an optical medium 
introduces additional possibilities. Indeed, changing the periodicity and material contrasts of 
the medium may permit tuning of the dispersive properties, e.g. the length scale on which 
a soliton can form may be altered. Thus, nonlinear and periodic structures are natural 
candidates for device design and applications. An example is the formation of centimeter- 
scale gap solitons in periodic optical fiber gratings. Such states have been shown to propagate 
at a fraction of the speed of light and have been proposed in schemes for optical storage and 
buffering; see, for example, 
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In the simplest setting, nonlinear electromagnetic waves in a one-dimensional periodic 
structure are governed by a nonlinear Maxwell equation: 

d 2 t (n\z)E + X \E\ 2 E) =d 2 z E. (1.1) 

Here, \ > is the Kerr nonlinearity coefficient, [3]. We assume a low-contrast, periodic 
refractive index profile, n(z), with mean n , given by 

n(z) = n + eN{z), n > 0, N(z) = N(z + 2vr), < e < 1; (1.2) 

n(z) is real-valued; no energy-dissipation has been included. The fluctuating part of the 
refractive index, N(z), can be expanded in the Fourier series 

N(z) = n p^ pz > N -p = N P , pe Z. (1.3) 



For simplicity, let us assume A^2 ^ 0. Then strong dispersion is excited by initial condi- 
tions of wave-packet type, i.e. a slowly modulated plane wave of a single frequency, chosen 
to be in (Bragg) resonance with the 7r-periodicity of the medium: 

E(z,t = 0) = e^ [E? (ez, 0)e iz + E{(ez, 0)e~ iz + c.c] , (1.4) 
where E 1 (Z,0) are spatially localized in Z = ez. This resonance strongly couples backward 



and forward propagating waves. In the choice of initial condition (1.4), dispersive effects 
which are set by the medium contrast, of size 0(e), have been balanced with nonlinear 
effects, by choosing the amplitude to be of size 0(e*). 

Suppose we make a formal multiple scale expansion based on the ansatz: 

E(z, t) = ea [E+(Z, T)e i{z - V ^ + E^{Z, T)e-' l{z+v ^ + c.c. + 0(e)] , (1.5) 
T = et, Z = ez, v g = l/n 
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Then if we only account for the principal harmonics, we shall arrive at the nonlinear coupled 
mode equations (NLCME) for Ef(Z,T): 



d T E? + v g d z E+ = iv 2 (N E+ + N 2 E^) + iT (\E+\ 2 + 2\E^\ 2 ) Ef 
d T E- - v g d z E^ = iv 2 g (N 2 E+ + N E~) + iV (\E^\ 2 + 2|£+| 2 ) 



(1.6a) 

;i.6b) 



where V = 3x/(2nl). Ef denote slowly varying forward and backward wave amplitudes; 
see [6] and references cited therein for details. 

NLCME has been rigorously derived as a leading order model in numerous contexts. For 
one-dimensional propagation of electromagnetic waves in nonlinear and periodic media, it 
was rigorously derived from the anharmonic Maxwell-Lorenz model in 12 . Derivations from 



the Klein- Fock as well as Gross-Pitaevskii equations have also been obtained; see 13 14 16 
17 . Explicit localized stationary solutions, called gap solitons, for NLCME are given in 
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The linear stability of the gap solitons was studied in p), and a linear, multi-dimensional, 
analog of NLCME was examined in j2]. 

However, NLCME is not the correct mathematical description of weakly nonlinear and 



weakly dispersive waves in the nonlinear and periodic Maxwell equation (1.1), (1.2). The 
deficiency of the NLCME system, (1.6), stems from the unperturbed (e = 0) equation being 
the non- dispersive one-dimensional wave equation. Due to nonlinearity, a single frequency 
initial condition, (1.4), excites infinitely many resonances, since e im ( z±t / n o) ^ m g % a \\ \[ e 



in the kernel of the unperturbed operator, n^d 2 — d 2 . In contrast, other models, such as 
the aforementioned anharmonic Maxwell-Lorenz system and the Gross-Pitaevskii equation, 
remain dispersive in the e = limit; this precludes infinitely many resonant modes. 

In [19], nonlocal equations derived from nonlinear geometrical optics and an equivalent 
system of infinitely many coupled PDEs, which take into account the infinitely many reso- 
nances, were systematically studied. One begins with the general weakly nonlinear ansatz, 



E{z, t) = e 1 * [E + (Z, T,z- v g t) + E~(Z, T,z + v g t) + 0{e 



which need not be nearly monochromatic. A necessary condition for the error term in ( |1.7 ) 
to be of order e on the time interval < t < O (e _1 ) is that the forward and backward wave 

) ±) ) 0± = z T v g t, satisfy the system of nonlocal evolution equations: 



components, E ± (Z, T, > 



{d T + v g d z + v 2 N^)E~ 



(d T - Vgdz ~ v 2 N d^)E p 



v g dj> 



2ty 



N(<p + 6)E-(Z, T,(f) + 29)de 



2tt 



2tt 



\E~(Z,T, 



de £ H 



N((f) - 6)E + (Z,T,(j) - 2d)dd 



(E 



-^3 



2ty 



\E+(Z,T, 



d9 E~ 



1.8a) 



;i.8b) 
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While we have omitted the ± subscripts on cL derivatives for the sake of brevity, the reader 



should note that in recovering the primitive field, as in (1.7), E + must be evaluated at <p + 
and E~ must be evaluated at </>_. E ± (Z,T,(p±) are assumed to be 27r-periodic in their <ft± 
arguments. A similar, but more general system of integro-differential equations was obtained 
though in that work, the authors set N 



m 
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and v„ 



If we expand E (Z, T, (j>) in a Fourier series with respect to the phase variable 



E ± (Z,T,<p) = J2E^(Z,T)e 



D ±ip<p 



1.9) 



the nonlocal system (1.8) may be re-expressed as a system of infinitely many nonlinear 
coupled mode (differential) equations for the Fourier mode coefficients, indexed by p G Z: 



d T E+ + v g d z E+ = ipv 2 g (N E+ + N 2p E~) 



+ iv— 
^3 



Y,E+EtE+ +r _ p + 2,[Y)E, 



-|2 
Q 



E: 



.q,r&Z 



3 T E~ - v g d z E~ = ipv 2 g (N_ 2p E+ + NqE;) 



+ iP 3 



^2 E q E r E q+r 



.q,r& 



(1.10a) 



(1.10b) 



In 19 the infinite system of PDEs (1.10) is referred to as the extended nonlinear coupled 



mode equations or xNLCME. Thus xNLCME is an extension of the classical NLCME (1.6) 



appropriate for highly resonant settings, such as the weakly periodic and nonlinear Maxwell 
model (O). Truncation of xNLCME to a single mode, Ef(Z,T), yields NLCME, (O), 



which, as noted, has spatially localized gap-soliton solutions. 

Numerical simulations of the primitive nonlinear and periodic Maxwell's equations, ( 1.1 ), 
give evidence of two phenomena. First, there appear to be long-lived spatially localized 
coherent structures. Second, within such spatially localized structures, a train of carrier 
shocks can form. These structures appear to be well described by xNLCME, 
The nonlinear Maxwell equation, 



(1.1) 
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dispersion which, as in the anharmonic Maxwell-Lorentz model 12 



does not incorporate any effects of chromatic 

takes off resonance the 



higher harmonics. However, chromatic dispersion on the length scales of many experiments 
is a negligible effect, j9J. Moreover, there are experimentally realizable regimes in which 
pulses with spectral content near the zero dispersion point are propagated [l5j. In these 
experiments, a broad band super continuum is generated. The carrier shocking mentioned 
above is a possible source of such broad band emission. 

In this paper, we explore, by analytical, asymptotic and numerical methods, the existence 
and properties of spatially localized structures of xNLCME. These coherent solutions have 
a full spectrum of active temporal frequencies and we therefore refer to them as broad band 
solitons. An earlier step in this direction was taken in 20 , where the authors studied 
what amounts to a truncation of xNLCME to first and third harmonics. Studying the 
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problem numerically, they found evidence for spatially localized solutions that they called 
polychromatic solitons. 



We focus on the stationary, small amplitude, near band edge, approximation of xNL- 
CME for a particular refractive index consisting of an infinite periodic array of Dirac delta 
functions. In this regime, xNLCME is well-approximated by an infinite system of coupled 
nonlinear Schrodinger equations, the extended nonlinear Schrddinger system, xNLS. We em- 
bed xNLS in a one-parameter family of equations, xNLS e , which continuously interpolates 
between a system of infinitely many decoupled NLS equations (e = 0) and xNLS (e = 1). 
Using bifurcation methods, based on the Lyapunov-Schmidt method and the implicit func- 
tion theorem, we prove the existence of solutions for a range of e G (— eo, eo)- By numerical 
continuation method, we establish the persistence of certain branches all the way to e = 1 
for finite truncations of xNLS e . Finally, we perform time- dependent simulations of xNLCME 
and find the small amplitude near band edge gap solitons to be robust. 

Outline of the paper: In Section [2j we present a direct derivation of xNLCME in the 
case of a periodic medium and show the sense in which xNLCME is an infinite dimensional 
Hamiltonian system. In Section [3] we heuristically determine conditions on N(z) for which 
we may expect exponentially localized gap solitons. This motivates us to focus on the case 
where N(z) is a periodic array of delta functions. 

In the small amplitude, near band edge, limit where xNLCME reduces to xNLS, we 
conjecture that localized stationary solutions of xNLS exist. Subject to this assumption, we 
prove in Theorem [T] that the gap soliton persists within xNLCME, in the asymptotic limit. 
Since the energy of xNLS is bounded below, it is natural to ask where a ground state of 
xNLS can be constructed variationally. Unfortunately, standard methods to not apply due 
to a loss of compactness, illustrated in Section 3.2.3| The existence of nontrivial critical 
points is an open problem. 

We therefore seek to construct localized states via a continuation method. First, we 
embed xNLS in a one-parameter family of systems, xNLS e , with e = corresponding to an 
infinite system of decoupled NLS equations and e = 1 corresponding to xNLS, the system of 
interest. In Theorem [2j we prove the existence of gap solitons for xNLS e for an open interval 
of |e| < e about e = 0. 

We next attempt to numerically continue xNLS e solitons on the interval [0,1]. In or- 
der to implement the numerical continuation, we seek approximate critical points of the 
xNLS e variational problem. To motivate this, in Section |4j we replace the variational char- 
acterization of xNLS 6 solitons by a finite dimensional minimization problem over families of 
Gaussian trial functions. We find critical points, with sign alternating amplitudes, of such 
finite dimensional approximations and give convincing numerical evidence that some can be 
continued to e = 1. 

In Section [5] we compute soliton solutions of truncated xNLS using information gleaned 
from the trial function approximations, and show that they are robust in time-dependent 
simulations of truncated. Section [6] summarizes our findings and highlights open problems. 
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2 Coupled Mode Equations 



In Section |2.1[ we present a derivation of xNLCME from Maxwell's equations using Fourier 
expansions of E ± (Z,T,(j)±), in the case where E ± (Z,T,<fr±) are periodic in <p±. In Section 



2.2[ we demonstrate that xNLCME is an infinite dimensional Hamiltonian system with two 
conserved quantities. 

2.1 Derivation of xNLCME in a Periodically Varying Medium 

For simplicity and without loss of generality, we set uq = 1 so that v g = 1. We rewrite the 



nonlinear Maxwell equation (1.1) with refractive index (1.2) as 

d 2 E - d 2 t E = 2eN(z)d 2 E + e 2 N(z) 2 d 2 E + X d 2 \E\ 2 E, 



(2.1) 



For e = and x = 0, (|2.1|) simplifies to the one dimensional wave equation with a solution, 

(2.2) 



given by the arbitrary superposition of right and left traveling waves, 

E {0 \z,t) =E + (z-t) + E-(z + t). 
For e small, we seek E = E e (z, t) in the form of a multiple scale expansion 
E{z,t) = e* (£ (0) {Z, T; z, t) + e£ (1) {Z, T; z, t) + 0(e 2 )) , 



(2.3) 



where Z = ez and T = et are slow spatial and temporal scales. Substituting (2.3) into (2.1), 
we obtain at first order in e: 

(d 2 - d 2 ) E^ = 2 (d t d T - d z d z ) + 2N(z)d 2 E^ + X d 2 \E^\ 2 E^. (2.4) 



The right-hand-side of (2.4) generates resonant terms along the characteristics of the wave 



equation, leading to secular growth of the correction E^ in (z,t). The slow evolution in 
(Z, T) is determined to remove these secular terms. 
We begin by expanding E in a Fourier series: 

E ± (Z, T; z,t) = J2 E^(Z, T)e ±ip ^, E^(Z, T; z,t) = £ Ej»{Z, T; t)e ipz . (2.5) 
Since E are real-valued, 



E±(Z,T) = E±(Z,T), peZ. 



(2.6) 
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Substituting (2.5) into (2.4), the terms of the equation proportional to e ipz are: 

(d 2 + p 2 ) E p = 2ip(d T + d z )E+e- ipt - 2ip(d T - d z )Ez/ pt 

+ 2J2<1 2 {N p - q E+e^ + N p+q E-e-^) 
i 

+ xJ2P 2 K E r^ +r - p e- ipt + 2 X £> ~ 2 <? + 2r) 2 E+E+E;_ r _ p e i ^+^ t 

q,r q,r 

+ x J2(p + 2q + 2r) 2 E-E;E + _ p _ q _ r e-^ +2 ^ t 

q,r 

+ x J2( P -2q- 2r) 2 E+E+E^ r e^- 2 ^ 

q,r 

+ 2 X J> + 2g - 2r) 2 E-E;E; +q _ r e-^- 2 ^ + x ^ E i E r E - P -i + r^ P \ 

q,r q,r 

where all sums are taken over Z. Removing the terms resonant with e ipt , we obtain 

(d T + d z )E+ =ip {n e; + n 2p e;) 



E+E+E*^ + 2E E t E t- 



q,r 



(2.7a) 



+ E E i E -i E -p + E £ K E U + 2 XXl 9 ^ 
i i i 

Removing terms resonant with e~ ipt , we obtain 
-{d T - 8 Z )EZ P =ip {N 2p EZ p + N EZ p ) 



+ ip~. 



r 



^2 E q E r E q+r+p + 2E E E « E i 



P+1 



q.r 



(2.7b) 



+ E t E -i E v + E o E E i E - P -i + 2 Y,\ E i\ 2E: 



where we have set T = 3 X /2 to be consistent with previous work, 12 19 . Exchanging p for 
— p in (2.7b), we have 

(d T - d z )E~ =i P {n_ 2p e; + n e;) 
r 



+ i p~ 



E E 1 E r E 1+r~p + 2E E E 1 E ~ 



P+1 



L 1,r 



+ K E -i E -v + E o E E * E v-i + 2 Y} E t\ 2E v 



At p = 0, (2.7) can be satisfied by choosing arbitrary functions Eq = Eq(Z =p T). For 
simplicity, we set Eq(Z =p T) = 0. If we additionally invoke complex conjugate relationship 
(2.6), (2.7) simplify to xNLCME, (1.10), from the introduction, provided v g = 1. 



Finally, the nonlocal system (1.8) can be recovered by introducing the identities 



E d 



(Z,T,<P) = ^^(Z.Tje^, E±(Z,T) = ±- f E ± (Z,T, 
P ez 71 J-7T 



(2- 



Constraints (2.6) imply that E ± are real- valued. Note that in the context of the primitive 
electric field variables, E ± (Z,T, (f)±) must be evaluated at different phases, <fr± = z t. 



2.2 Hamiltonian Structure of xNLCME 

Let Eq , and define H = J R T-idZ , where 

n = \ E - - E ;d z E p - e;o z e; + E-d z E-) 

1 p p 

- N ^T(|£+| 2 + \E;\ 2 ) - N 2p (El p E± p + E;E+) 



r 

6 



z2 E P E - P ) [zZ e p e - 



r 

6 



]• 21" 

~~ fi ( E p E q E r E t+r-p + E p E q E r E q+r-p) ~ IT 



Y) E l 



(2.9) 



-i+|2 



p,q,r 



with all sums are over Z \ {0}. Then, xNLCME has the structure of an infinite dimensional 
Hamiltonian system: 



in 6H 3 F~ 
' ip 5E+> 9tE ? 



-i P j^, P ez\{0}. 



(2.10) 



Formally, the Hamiltonian (|2.9|) is conserved under the flow of xNLCME. Besides the 

(2.11) 



Hamiltonian, the total power N = j R MdZ is invariant, where the density is 



pez 

This follows by direct computation. 

Since N 2p = N_2p, G Z, the symmetry of equations (1.10) implies that if the constraint 
Ep = E^ p , associated with real initial conditions for E , is satisfied at T = 0, then it is 
satisfied for all T. Additionally, if Er are zero initially for even p, they remain zero for all 
time. This allows us to restrict (1.10) to the odd harmonics, p G Z 0( jd, and set 



Et = 0, pez c 



(2.12) 



Under these constraints, the conserved integral (2.9) reduces to the Hamiltonian: 
U = \ E \{E;d z E; - E p d z E7 p - E;d z E; + E;d z E p ) 



2 ^ p 

pez odd 



-n E (i^pT + i^i 2 )- 2 E 

pez odd pez odd 



vpez odd 



(2.13) 



-|2 
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E ( E P E t E r E t+r-V + E P E 1 E r E q+r-p) 



p,q,r££ odd 



As in the case of standard NLCME, (2.13) is unbounded from above and below subject to 
the constraint of fixed M . Thus, critical points are expected to be of infinite index. This 
suggests that variational methods will be of limited applicability for studying the stability 
of localized stationary states of xNLCME. 



3 Gap Solitons 



We now begin to explore the existence of localized stationary states of xNLCME (1.10), 
called gap solitons. Setting v g — 1 for convenience, we seek solutions of the form 

E+{Z,T) = e i ^ N °-^ T A p (Z), E;{Z,T) = e^ N °-^ T B p (Z), p G Z, (3.1) 

where fl is a real frequency parameter and {A P (Z), B p {Z)} p ^i are complex-valued ampli- 



tudes. Using constraints (2.6) and (2.12), we assume 



A p — A_ p , B p — B_ p , p G Z dd, A p — B p — 0, p G 7L C 



(3.2) 



The infinite family of amplitudes {A p , B p } pe z odd satisfies the extended system of stationary 
equations 

iA'p(Z) + pQAp + pN 2p Bp + p^ hA p J2 \ B i\ 2 + E 4AA>-<z-r J = 0, (3.3a) 
-iB' p (Z) + ptlBp + pN 2p A p + p^ I 3B P ^ | A q \ 2 + £ B q B r B p _ q _ r J = 0, (3.3b) 



with constraints (3.2). Linearizing about the zero solution yields decoupled systems of dif- 



ferential equations with solutions 

A p 
B n 



,±Zy/\pN p \*-(pSl)* 



pez 



odd- 



(3.4) 



A sufficient condition for spatial localization near the zero solution is only possible if \Q\ < 
Qq = min pe z od al-^2p|, implying three possibilities: 
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Case 1, fi > ; An example would be N 2p = 1, p 6 Z, in which case the refractive index, 
N(z), is a periodic sequence of Dirac delta-functions. 

Case 2, Q = and min pgZodd \pN 2p \ > 0: An example would be N 2p = p^ 1 , p G Z odd , for 
which N(z) would correspond to a periodic sequence of step functions. 

Case 3, Q = and min pgZodd \p 2 N 2p \ < oo: In this case, N(z) is continuous. 

If N 2p = 1, p G Z odd , the band gap of each mode is opened, and the widths of the band 
gaps grow as \p\ — > oo. However, because of the coupling between the Fourier modes with 
amplitudes {A p , B p } pe z odd , the stationary localized mode (gap soliton) may only reside in 
the gap of a fixed width, < Q = 1. 

If N 2p = the band gap of each mode is again opened, but the widths are nearly 

constant as \p\ — > oo. However, the band gap for the coupled gap soliton shrinks now to zero 
and the parameter fl must be set to 0. 

If N 2p = <D(\p\- 2 ), the widths of the band gaps shrink with the larger values of p, and 



the exponential decay (3.4) ceases as \p\ — > oo, even if fi = 0. We do not anticipate the 
existence of gap solitons in this case. 

We now restrict our attention to Case 1: Q > and set N 2p = 1 for all p G Z odd . System 



(]3.3|) can now be rewritten as an equivalent integro-differential equation: 

1 



-d z + Sid*) A + d^B + -d 4 
(d z + m^)B + d^A + 



A d + 3 



B" 



2tt 
1 

2^ 



\B(Z, s)\ 2 dsj A 
\A(Z,s)\ 2 ds] B 



0. 
0. 



where we have introduced the two Fourier series, 

A(Z,<j>) = A p(Zy p ^ B(Z, 
pez dd 



1 P 



(3.5a) 
(3.5b) 

(3.6) 



We note that if one wishes to compute the primitive electric field induced by these envelopes, 
care must be taken in where the phase variable, 0, is evaluated. Indeed, the electric field 
associated with {A p , B p } pG % odd is given by 



E(z,t) 



e ip{N °- n)et e ip(z - t) A p (ez)+ ^ e ip{Na - n)et e-' ip{z+t) B p (ez) + 0(e) 



ip(N -n)et 

[A(ez, (N - tyet + z-t) + B(ez, (N - tyet - (z + t)) + 0(e)] 



(3.7) 



in agreement with the ansatz (1.7). 
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3.1 NLCME Gap Solitons 

As noted, the truncation of xNLCME to E 1 yields the classical NLCME. We now review 
the details of the NLCME gap soliton. 

The spatial profiles of NLCME's gap soliton are given by solutions of the stationary 
equations: 

iA[(Z) + QAi + B 1 + TdAI 2 + 2\B l \ 2 )A 1 = 0, (3.8a) 
-iB[(Z) + nB x + A X + T(2\A 1 \ 2 + \B 1 \ 2 )B 1 = 0. (3.8b) 

For Q E (—1, 1), these equations admit the exact solutions: 

= \Hacosh(^iZ)-if3smh(fIZ) , (3 ' 9a) 

Bi{Z) = 7" — (3.9b) 

V 3r a cosh(/iZ ) + lp sinn(/iZj 

where 

a = Vl + tt, /3 = y/l-n, fJi = Vl - tt 2 = a/3. 



The localized solution (3.9) satisfies the symmetry property 



A 1 (Z) = A 1 (-Z), Bx(Z) = B 1 (-Z), Ze 



We shall call the solution of (3.9) a monochromatic gap soliton, since the associated 
approximate solution of the nonlinear Maxwell model consists of a slowing varying and 
localized envelope with a single fast (carrier) frequency of oscillation. This is in contrast 
to the broad band, or polychromatic, solitons which possess slowly varying envelopes on 
multiple distinct carrier frequencies. It seems unlikely that there is an explicit solution of 



the system (3.3) of infinitely many coupled mode equations. 



3.2 Persistence of Solitons in a Band Edge Approximation 

We now explore a small amplitude, spectral band edge, approximation of xNLCME, which 
will lead to an infinite system of coupled NLS type equations, xNLS. 

3.2.1 The Band Edge Approximation 

The gap in the continuous spectrum exists for Q G (—1,1). In the truncated coupled mode 



equations (3.8), the exact solution (3.9) shows that the amplitude ||Ai||i,°o of the gap soliton 
becomes small as f2 — > 1. Using the parameterization f2 = a/1 — and the asymptotic 
expansion 

A l =^U 1 {Q + 0{f), B 1 = -//CMC) + 0(fi 2 ), 
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where ( = fj,Z is slow variable and \x is a small parameter, we can formally reduce the system 

(3.10) 



of differential equations (3.8) to the scalar second-order equation for t/i(C) : 

c/('(C) -CMC) + 61^(0 = 0. 



This equation admits the localized solution 



B*(C) 



1 



sech(C), 



(3.11) 



which corresponds to the asymptotic approximation of the gap soliton (3.9) as O — > 1. 



Generalizing this approach to the system of infinitely many coupled mode equations, 

(3.12) 



(3.3), we substitute = a/1 — /i 2 and 

A p = nAp(C), B p = -fiB p (C), p e Z odd 



with C, = into the coupled mode system (3.3) to obtain 



ifiA' p + pV 1 - ^ A p - P B p + fr/i 2 F p = 
ifi&p - p^l - y?B ¥ + pA p - \Yn 2 G p = 0, 



where 



f p = 3A P J2 + Yl K~ A A 

G p = 3B p \Aq\ 2 + 2^ B q B r B 
<?6Z odd 

Introducing the variables 



q,r£Z odd 



p—q—ri 



p—q—ri 



u„ 



A v + B v 



i r p 



A p — Bp 



the system (3.13) can be written as 



2pV p + ifiU p + ( Vl - /i 2 - lj pV p + -Tfi 2 p(F P ~ G p ) = 0, 
W'p + v 1 ~ ^ - 1 pU p + l -Tw(F p + Gp) = 0, 

where F p , G p are rewritten after the substitution of the new variables. 



Now, if we formally expand in powers of /i, 



U p = U p + 0(n 1 ), Vp = -^-U p + 0(n 2 ), 



(3.13a) 
(3.13b) 

(3.14a) 
(3.14b) 



(3.15a) 
(3.15b) 



(3.16) 
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we find obtain, at leading order, an infinite system of coupled NLS type equations, that we 
deem xNLS: 

U'^Q-tfUv+^-vLu, Kl 2 + E E W^-,-r)=0, pGZ. (3.17) 



d <?ez odd rSZ odd 

This can be rewritten as the integro-differential equation 



(d 2 c + dl)U= -Tdl 



after introducing the Fourier relations 

£/(C,0) = E ^p(Oe W , ^p(C) 



2^ 



|[/(C,0)| 2 d0 £/ 



(3.18) 



2tt 



C/(C,0)e- 



We will now justify the reduction to xNLS, (3.17) 



3.2.2 Preliminaries 

We first introduce appropriate function spaces in which we study the problem. Let T denote 
the interval [0, 2tt], with endpoints identified so that functions on T are understood to be 2n- 
periodic. We shall consider functions defined on R x T, admitting the Fourier representation 

f/(C,0) = E f/ P (Oe i ^. 
For any s, the function space X s is defined by 



X s = jf/(C,0) e H S (R x T) 
and equipped with the norm 

\\U\\xs = 



u(c,<f>) = u({,<f>), 

U(C, 0) cos(2 P 0)c/0 = 0, VC e E, p e N 

\ 1/2 



(3.19) 



E f(p 2 +er\u P (o\ 2 dA 



(3.20) 



We shall frequently go back and forth between the U and {U p } pe z odd representations of 
functions in X s . 

The Sobolev space H S (WL x T) is a Banach algebra with respect to the pointwise multipli- 
cation for any s > 1. Moreover, from the continuous embeddings H S (M> xT)4 L°°(IR x T) 
for s > 1 and Z 2 (Z) e — >■ Z°°(Z), we infer that if U G X s for s > 1, then 



lim t/(C, 



0. 



e T. 



(3.21) 
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Let Bs(X s ) denote a ball of radius 6 in Banach space X s centered at the origin. The 



Hamiltonian H with the density (2.13) consists of the terms controlled by the H 1 norms of 
E ± . To see this, recall the continuous embedding if 1 (IR xT)4 L A (R x T). It follows that 
for any E ± G Bs(X s ) with s > 1, there is a constant Cg, s > such that 



H< C* 



X' 



+ \\E- 



x 



Furthermore, the map (E + ,E ) h- >■ H is continuous in X s . Although we will mainly study 
the problem in X s with s > 1, we note that the energy is well defined in X . 



3.2.3 Proof of Result 



We now rigorously justify the small amplitude approximation of (3.3) by (3.17) 



Theorem 1. Fix s > 1 and assume the existence of localized solution U G X s to (3.18) 
satisfying the reversibility symmetry, 



U p (() = U p (-(), p6Z odd , Ce 



(3.22) 



Also assume that the linearized operator of system (3.18) at U is invertible in the subspace 



of X s associated with the constraint (3.22) 



There exists /i > such that for any fi G (—fio, Ho)> the system of stationary coupled 
mode equations (3.3) with = a/1 — admits a unique localized solution A,B e X s 
satisfying the symmetries, 



and the bound, 



A P {Z) = A p (-Z), 



\A — \iU {\x- 



B p {z 



B p (-Z) } p G Z odd , Z G 

X s + \\B + fiU(fi;-)\\ X s <CfI 2 . 



(3.23) 
(3.24) 



Proof. First, we note that the vectors fields F(A,B) and G(A,B), defined by their compo- 
nents in (3.14), are analytic (cubic) maps from X s x X s to X s for any s > 1. Eliminating 



U p from system (3.15), we obtain 



v 2 % - v P " = - 6 r 



v 



^-i)(f p -g p )- w (f' + g;; 



any s > 1, where the X s 



(3.25) 



for 



The right-hand side of system (3.25) defines an analytic (cubic) map from X s to X s 

" 2 norm is of order O(fi) as fi — > oo. The left-hand side operator of 

' 2 to X s 



system (3.25) has a bounded inverse from X s 
in X s 



thanks to the zero mean constraint 
By the Implicit Function Theorem, we infer that for any 5 > and any s > 1, there 
is /i > such that for all /i G (— /i ,yUo) an d for all U in a ball Bs(X s ), there is a smooth 
map X s 3 U i— > V[U] G X s which solves system (3.25) and satisfies the bound, 



3C > : 



\V\ 



\x> 



<Cfi, fie (-/lo, Ho), U G Bs{X' 



(3.26) 
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On the other hand, eliminating V from system (3.15), we obtain 



U" 

v 



p 2 u p + f 



^ + l)(F p + G p )-tfip(F;-G' p 



0. 



(3.27) 



Thanks to the bound (3.26), the cubic terms of the system (3.27) are different from those 
of the system (3.17) by the error of the order of 0(fi 2 ) in X s-2 . Under the assumptions of 



the existence of the solution U £ X s of the truncated coupled NLS equations (3.18) and the 



invertibility of the linearized operator in the subspace of X s associated with the constraint 
(3.22), the linearized operator has a bounded inverse from X s ' 2 to X s for any small 



By the contraction mapping arguments, there is a solution U near U in X s such that 

3C>0: \\U- U\\ x > < C/j 2 . 



This gives the statement of the theorem, after the original variables A, B, and Z are restored 
from the transformations above. □ 



3.3 Hamiltonian &; Power of xNLS 



The extended system of coupled nonlinear Schrodinger equations (xNLS) (3.17) inherits the 
Hamiltonian structure of the coupled mode equations (3.3). The energy functional for (3.17) 
is given by 



xNLS 



1 



.77' z -I- 77 

2 1 u p | -t- \u p \ 



_pez odd 

We also define the power, 



vpez odd 



r 

3 



p,9,reZ odd 



X 



xNLS 



.pGZ 



d(. 



d(. 

(3.28) 
(3.29) 



Energy functionals are often used in proving the existence of localized solutions to con- 
strained variational problems, e.g. 



minimize i/xNLS subject to fixed iV X NLS- 



(3.30) 



Unfortunately, this strategy fails for our problem, as demonstrated by the following coun- 
terexample. Let 

U P (C) = K/ 2 W(XnC) (S p ,n + Sp,-n) , P G Z odd , (3.31) 

where W G i/ 1 (lR) is a fixed function, A n > is an arbitrary parameter, and n > 1 is an 
arbitrary odd integer. Then, iV xN Ls is independent on the parameters A n and n. On the 
other hand, 

2A 2 

H x nls = — Wl 2 ~ 6A n ||£?|| L 4. 



If we set X n = n and let n — > oo, we obtain no lower bound on if x NLS- Thus, localized 
solutions of xNLS, (3.17), if they exist in some X s , cannot be global minimizers; they will 
either be local extrema or saddle points. 
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3.4 Persistence of Monochromatic Solitons to Coupling in xNLS 

We study the question of persistence of NLS solitons within xNLS by embedding xNLS in a 
one-parameter family of models, xNLS e , for which xNLS is an infinite system of decoupled 
NLS equations and xNLS 1= xNLS. Our formulation is: 



U p \() - p 2 U p + 6p 2 TU^ + ^eT [3U P £ \U q \ 2 + ]T U q U r U t 



p-q-r 



0, p G Z, 



odd; 



q,rE& odd 



(3.32) 

where the sums exclude the cubic self interaction terms, Up. Within each mode of the 
decoupled system at e = 0, we have a solution 

U p {() = ±C/*(pC), P e Z odd , (3.33) 



where U*(() is the NLS soliton ( |3.11p 

We now prove the persistence of (3.33) within xNLS e , (3.32), for all e sufficiently small. 
Without loss of generality, we can take p = 1. Furthermore, we make the reduction 

U p (() = U P (C) = U_ p ((), , P eZ oAd . 

In other words, we now assume that the envelopes in each harmonic are real-valued. 

Theorem 2. Fix s > 1. There exists eo > and C > such that for any e G (— eo,eo); 
xNLS, (3.32) admits a unique localized solution U G X s satisfying the even symmetry: 

U P (C) = U p (-C) (3.34) 
Moreover, U ((, <fi) is a small deformation of the unperturbed e = soliton solution: 

Z4(C,0) = 2tf*(C)cos(0), 

in the sense that 



\U-U\ 



X' 



< Ce 



(3.35) 



Proof. The proof relies on a Lyapunov-Schmidt reduction where we shall first express the 
higher harmonics as functions of U% = C/_i, and then apply the implicit function theorem to 
an equation written entirely in terms of U\. 

From (3.32), define F in terms of the components 

F p = 3U p l^! 2 + U « U r U i 



p—q—ri 



where each F p excludes the purely self-interacting terms. For |p| > 1, we can clearly write 



U P = ~{dl - P 2 )-y (6rt/ p 3 ) - e(<9 ( 2 - p 2 )-y 2 fF p . 
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(3.36) 



The terms on the left are in X s since s > 1 and (d 2 — p 2 )~ 1 p 2 is a bounded operator. 
Therefore, for sufficiently small eo > and finite Sq, the contraction mapping theorem yields 
a unique map, 

$:([/ 1 , e )^{C/ p } |p|>1 

in a ball U\ G Bs(H s ) with 5 < 5 and |e| < e . For a given f/i = LLi, we have expressed 
the other modes in terms of this fixed profile. Form (3.36), we can see that there exists a 
constant C > such that for all |e| < eo, 

H^t/i^lU^CellC/Jip 



We now eliminate {f^}| P |>i from the p = ±1 equations of (3.32) using the mapping $. 
Since U\ = Z7_i, we only consider the p = 1 equation: 



u" - u x + ert/ 3 



;erF 1 [$(f/ 1 ,e)]. 



(3.37) 



For any U\ G B$(H S ) with finite 5 > and small |e| < eo, there is a C > such that 

\\F x [®(U u e)]\\ B .<Ce\\Uif Ht . 



To solve (3.37), we hope to apply the implicit function theorem. Thus, we must linearize 



(3.37) near U* at e = 0, and show that it is invertible. The kernel of the linearized operator, 

c\ 2 - 1 + 



is just dcJJ*. Obviously, this does not satisfy the symmetry constraint, (3.34). Subject to 
this condition, the operator is an isomorphism from H^ ven — > H^~^, the subspace of H s of 
even functions. Hence, the implicit function theorem yields a neighborhood of U+ in H s in 
which we can obtain U\ with |e| < e x < e . 

Moreover, we see from (3.37) that there exists C > such that for all |e| < e 1; 

\\Ui-UJ\h° <Ce 2 



Combining this estimate with (3.36), yields (3.35). 



□ 



This result has several obvious extensions. We can consider the local continuation about 
a soliton localized in any other mode po G Z dd, 



u P0 (() = ±uM)- 



(3.38) 



We could also continue a solution about any finite collection of such solitons. However, if we 
begin with solitons in every odd harmonic, they will not have finite L 2 , as 



/ £ \u P \ 2 d( = \\u,\\ 2 L , £ i 

pez odd p6Z odd 



-...i-i pez odd ^ 

which diverges. The continuation of such infinite energy solutions is an open problem. 
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4 Variational Approximations 



As noted at the end of Section 3.2.3 although the functional is bounded from below, the 
natural variational formulation for localized solutions of xNLS exhibits a loss of compactness. 
In this section we explore the use of this functional to obtain Rayleigh-Ritz or Galerkin-type 
approximations to such states. The parameters in these approximations can be uniquely 
determined from the conditions that the ansatz gives a stationary point of i? x NLS- Indeed, 
variations of -£f x NLS produce the Euler-Lagrange equations, which are equivalent to the dif- 
ferential equations (3.17). 



4.1 Gaussian Approximations 

Let us consider the Gaussian variational ansatz 

U P (0 = a p e- b ^'\ peZ odd , 



(4.1) 



where a p G K and b p G M. + are parameters of the variational approximation. The Gaussian 
ansatz is useful because all integrals in i? x NLS can be computed in the analytical form. Direct 

substitution and integration show that J --HxNLS becomes 



H, 



Gauss 



E 



a: 



b a 2 



r E 



a p a q 



¥ E 



y/2, 



a p 



p,g,re^, odd 



yjb p + b q + b r + b p _ q _ r 

(4.2) 

If the artificial small parameter e is introduced to decouple the different modes, as in 



system (3.32), then (4.2) is rewritten in the form, 

~2 



H, 



Gauss (^ J — 



E 



b p a 2 p 



P 



-*[ E 



a P a q 



3at 

r p 



I E 



V2a p 



(4.3) 



p.gt^odd 



V b P + b 1 3 M ,^t dd V b P + K + K + b p - q - T 



The above sums with e as a prefactor exclude the purely self interacting a*/ ' \fb~ p terms. 

If e = 0, there exists an uncoupled solution of the Euler-Lagrange equations produced 
from variations of ifcauss(O), 



2 3/4 

: 3fV2' 



, pel 



odd- 



(4.4) 



The exact solution (4.4) will be used as a seed point in the numerical continuation algorithm. 
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Table 1: Computed values for a truncated Rayleigh-Ritz approximation for e = 1. 



No. of Modes 


a\ 


6i 


a 3 


bs 


a 5 


h 


1 


0.56060 


0.33333 










2 


0.56321 


0.33148 


-0.13918 


3.9413 






3 


0.56329 


0.33189 


-0.14585 


3.6287 


0.062822 


8.5577 



4.2 Numerical Continuation 



Truncating H Ganss (e) in (4.3) to resolve only iV harmonics, we define i^Q auss (a, b, e). The 



associated system of 2N Euler-Lagrange equations is 

V a ^ auss (a, b, e) = 0, V b ^ auss (a, b, e) = 0. (4.5) 
We now seek solutions of the e = system, where all modes are decoupled, that can be 



continued to e = 1, the desired system. The natural family of solutions is given by (4.4). 
Thus, for our e = starting point, we consider solutions of the form 

2 3/4 

a * = jjfya Oi^-.-iO-jv) (4.6) 

b, = i(l 2 ,3 2 ,...,(2iV-l) 2 ) (4.7) 



where o~j G { — 1,0,1}. The variances, b*, are unaffected by cr. Indeed, for cxj = 0, b 3 



is 



ill-defined, and can take any value. 

We now explore continuations from various cr's. Before giving the results, we state the 
conjecture that our computations suggest: 

Conjecture 4.1. For any N > 1, there is a nontrivial configuration cr that can be continued 
from e = to e = 1. At e = 1, the amplitudes are sign alternating, 



si 



ign(a p ) = (-l)<W-U/2. 



For a system of two modes (N = 2), the numerical continuation of four cr configurations 
is plotted in Figure [TJ The configurations cr = (0, 1) and cr = (1, —1), can be continued to 
e = 1, while the other two collide and terminate near e = 0.368. Extending this to a system 
of three modes, we plot the analogous results in Figure |2} Three configurations cr = (0, 1, 0), 
cr = (0, 0, 1), and cr = (1, —1, 1) can be continued to e = 1. We note that the configurations 
cr = (0, 1), cr = (0, 1,0), and cr = (0,0, 1) are trivial in the sense that they are generated 
by the reductions of the truncated coupled NLS equations. When more modes are included 
into the system, these degenerate configurations are destroyed. On the other hand, the 
configurations cr = (1, —1), (1, —1, 1) are non-trivial and persist with respect to adding more 
modes in the coupled NLS equations. Our results for the non-trivial configurations at e = 1 
are summarized in Table [U 

Though computations for two and three modes suggest that an alternating configuration 
of ±l's can always be successfully continued to e = 1, this is not the the following 
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0.8r 
0.75 

0.7- 
0.65- 

0.6 
" 0.55 - 

0.5 
0.45- 

0.4 
0.35 



^ , - r-> , - : ^ 



... (1.0,0) 
-(1,1,0) 
(1-1,0) 

-(1,0,1) 
(0,1,1) 
(0,1,-1) 

-(1,1,1) 
-'(1,-1,1) 
' (1,-1,-1) 

-(0.0,1) 
-(0.1,0) ■ 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



1 




.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 




Figure 2: Various continuation branches for a three-mode Rayleigh-Ritz system. 
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Table 2: Computed values for a truncated Rayleigh-Ritz approximation with fixed b* 
for e = 1. The branch from which we continue is alternating for 1 < iV < 4: 
cr = (1), (1, — 1), (1, — 1, 1), (1, — 1, 1 — 1). The case iV = 5 is continued from the branch 
(1, —1, 1, 1, 1) and N = 6 is continued from (1, —1, 1, —1, —1, —1). 



No. of Modes 


a % 


a 3 


a 5 


a 7 


a 9 


an 


1 


0.56060 












2 


0.5643 


-0.12734 










3 


0.56409 


-0.13759 


0.061454 








4 


0.56386 


-0.14037 


0.068618 


-0.037695 






5 


0.56372 


-0.14139 


0.071254 


-0.042822 


0.026041 




6 


0.56364 


-0.14184 


0.072457 


-0.045015 


0.029896 


-0.019323 



computations demonstrate. We first make the additional simplification, observing that the 
values of bj in Table [l] are close to j 2 /3. This motivates fixing them as such, and only solving 
the problem for the amplitudes, a. Thus, we solve 



V.fl£ n Ja,h.,e) = (4i 



where b* is given by (4.7). The results of our computations with these fixed variances are 
given in Table [2] Continuation from the alternating branch cr = (1, — 1), (1, — 1, 1), . . . is 
successful till N = A. The alternating branch cannot be continued to e = 1 for five and six 
modes, though there are other initial states that can be continued to e = 1, with sign sign 
alternations at e = 1; see Table [2j 

Though these results were initially computed using a naive continuation algorithm in 
Matlab, solving with a given value of e and using that solution as the initial guess for a 
larger value of e, they were confirmed by our computations using AUTO [7j[8]. 

Though the starting branch may not have an alternating sign structure, sign alternating 
solutions may still be found at e = 1. This makes it challenging to perform numerical 
continuation with these branches if we no longer assume the variances to be fixed. For a 
system of five modes, aj must change sign. When it crosses zero, the variance becomes ill- 



defined introducing numerical difficulties. On the other hand, if we iterate the system (4.5) 



near the solution of (4.8) for e = 1, the convergence is usually achieved with few iterations. 



4.3 Tails of the Variational Solutions 

Though we are able to construct a sequence of Rayleigh-Ritz approximations with Gaussian 
ansatz, it is not yet clear if such solutions should exist in space X s for s > 1 or at least have 



finite power (L 2 ) in the limit iV — > oo. Indeed, the solution (a*, b*) given by (4.4) for e = 
with all a p ^ has infinite power, since 



\a p exp(-6„C 2 )|^C 



2 x3/2 



3J |p|P 
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and J2 P ez d d R = °°' However, the results of Table |2j show that at e = 1, the sign-alternating 
amplitudes {a P } p ez odd are also decaying in p e Z odd . We explore whether or not the decay 
is sufficiently rapid to have finite power and to belong to the energy space, where H G3bUSS is 
finite. To this end, we employ a more refined trial-function ansatz, allowing for weak decay 



of a p : 



a p = A(-l)^-^ 2 \p\-\ b p = P Y , P 6Z odd (4.9) 

where A and 7 are unknown parameters to be determined from the Euler-Lagrange equations. 
If 7 > 0, the Rayleigh-Ritz approximation has both i/causs and iVcauss finite. 



Substituting (4.9) into (4.2) yields a two parameter Hamiltonian 

h{ 1) A) = A 2 f{ 1 )-A i Yg{ 1 ), 



(4.10) 



where 



f(l) 
0(7) 



E 



\p\ 



-1-27 



+ 



E ^ 
E 

p,q,rEZ odd 



■sjp 2 + q 2 

_]_^([p[+|?[+|r[+|p-9-r|)/2 



|P|" 



7 |?| 1 r l J \P 



\Jp 2 + q 2 + r 2 + (p — q — r) 2 



Solving 9,4/1(7, A) = 0, we find 



Plugging back in, we get 



A 2 ( 7 ) 



2^(7) 



%) = / i ( 7 ,A( 7 )) 



f{if Sin? = 1 f{if 

2Tg( 1 ) 4^(7) AT g(>y) 



(4.12) 



Truncating this approximation to N modes, ft, (7), we are able to identify a sequence 
of critical points, suggesting convergence as N — > 00 and the existence of a critical point in 



the primitive functional, (4.12). A few of these approximations are plotted in Figure [3] with 
T = 1. All of the computed/?^ (7) 's have the property that 



lim ^(7) = ^(7) 



7— >oo 



9 V 3 



(4.13) 



The critical values of 7, 7*, are given in Table [3j These appear to converge to a value of 7 
near 7 = 1.26 indicating that the corresponding variational approximations belong to the 
energy space of the coupled NLS equations. Moreover, since 



pez odd 



2s-li 1 2 

\ a p\ ~ 



E ipi 

pez odd 



2s-l-27 
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Figure 3: Two-parameter approximation (4.9) of Hs auss for different truncations 



Table 3: Computed critical values of 7 for the curves in Figure [3J 



No. of Modes 


7* 


A 7 * 


2 


1.35511 




3 


1.30184 


0.05327 


4 


1.28176 


0.02008 


5 


1.27208 


0.00968 


6 


1.26672 


0.00536 



and 7 ~ 1.26, the corresponding variational approximations belong to the space X s for 
s < 7. Therefore, the results of Theorems [T] and [2] can be used in the nonempty interval for 
the values of s G (1,7). As it appears that 7 is strictly greater than one as the number of 
modes increases, a solution of infinitely many modes might be more regular than H l ; indeed, 
it would be Holder continuous. 



The sign alternating structure of the ansatz (4.9) is fundamental for the existence of the 



critical point of hfa, A). For the variational ansatz, 



A\p\- 



el 
3 ' 



pez, 



odd; 



(4.14) 



we can redo the computations to obtain Figure 111 No critical point of h("f, A) exists for the 



sign-definite variational approximation (4.14). 



5 Numerically Computed Gap Solitons 

Using our observations from the Rayleigh-Ritz approximation, we are motivated to solve the 



xNLS, (3.32), directly for existence of the gap solitons. We note that in 20 , the authors 



explored the related problem of broad band solitons of xNLCME truncated to two modes. 
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5.1 Computation of the Gap Solitons 



We numerically solve equations (3.32) by continuation. Our starting point is the exact 
solution at e = 

f/ p (C;e = 0) = ^sech(K), 



where er is a branch found in Section |4~2 that led to a non-trivial solution at e = 1. Iterating in 
e, we solve the system (3.32) using Matlab's bvp5c algorithm with absolute tolerance 10 -4 , 
relative tolerance 10" 

We use the even symmetry of the 



" 8 , on the domain [0, 25]. bvp5c is a nonlinear finite difference algorithm 
for two-point boundary-value problems discussed in 



18 



solutions to impose the boundary condition U'(0) = 0, and the artificial boundary condition 



^p(Cmax) +pUp{( n 



0. 



The results for systems of up to six coupled NLS equations at e = 1 appear in Figure [5j As 
we can see the amplitude decays in p, and they appear to approach some asymptotic profile. 
We conjecture that this profile persists as additional modes are included. Alternatively, the 
solution can be expressed as U((, 9) by combining the Fourier modes. The resulting solution 
surface of the integral-differential equation (3.18) appears in Figure 16} The inclusion of 



additional harmonics induces a more ornate structure near the extrema. 

Though we have computed these finite truncation solutions, we reiterate the question 
whether the corresponding solutions have finite power. For our computed solutions, we find 
that the power, A^nls, appears to converge and most of the power remains in the first mode. 
The data is given in Table |4} 



5.2 Gap Solitons in Time Dependent Simulations 



Small amplitude gap soliton solutions of the coupled NLS equation (3.17) can be used as ini- 
tial conditions in the coupled mode equations (1.10) to assess their stability and robustness. 
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Figure 5: Soliton profiles for the coupled NLS equations (3.17). 



Table 4: Computed powers for the soliton profiles appearing in Figure [HJ 



No. of Modes 


\\Ui\\h 


|A^ X NLS 


1-ZVxNLS ~ 


Ui\\h 


1 


0.66667 


0.66667 







2 


0.66982 


0.68582 


0.016000 




3 


0.67147 


0.68929 


0.017825 




4 


0.67211 


0.69031 


0.018201 




5 


0.67226 


0.69070 


0.018441 




6 


0.67236 


0.69088 


0.018523 
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-4 C 

(a) Two Modes 



-4 C 
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Figure 6: The solution surface of the integral- differential equation (3.18) generated by the 
truncated coupled NLS soliton on Figure [5j 



Once the solution {U p (()} p( zz od 
simulation are given by 

El 



is computed, the initial conditions for the time dependent 

. p \Z,Q) = nU p {p.Z), E;(Z,Q) = -nU p {p,Z), p e Z odd , (5.1) 

with even modes set to zero. By Theorem [TJ the small amplitude approximation is only 
accurate up to 0(fi 2 ). We explore this small error as a source of the initial perturbation. 
We present the results of two and four mode systems. In each case, we truncated both 



the the system of coupled NLS equations (3.17) and the coupled mode equations (1.10) at 
the same number of resolved modes. In our simulations, we take as our constants 



1, N = 0, N.. 



2p 



1. 



The simulations were performed with the indicated number of grid points using a pseudo- 
spectral discretization and RK4 time stepping. For both the two and four mode simulations, 



the initial conditions (5.1) were computed with greater precision than an in the previous 
section; the absolute tolerance was 10~ 7 and the relative tolerance was 10~ 9 , and the domain 
was [0, 35]. 



In Figure |7| we plot the normalized time-space surfaces of \E£\ from our simulations of 
the first four odd modes. For both values of /j,, the solution is persistent, but the oscillations 
are greater for the larger value of fi, and there is some decoherence near the peak. With 
the smaller value of /i, there is far less distortion. Additional details of the dynamics are 
available online in the following animations: 

Two Mode Truncation The following simulations were computed with 1024 grid points. 
The /i = .4 simulations were computed on the domain [—50, 50], the \i — .2 simulations 
were computed on the domain [—100, 100], and the // = .1 simulations were computed 
on the domain [-200,200]. 
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http : / /www .math . toronto . edu/simpson/f iles/media/broadband/model_13 .mp4 



http : //www . math . toronto . edu/simpson/f iles/media/broadband/mode3_ 13 . mp4 



Four Mode Truncation The following simulations were computed with 2048 grid points. 
The /i = .4 simulations were computed on the domain [—50, 50], the /i = .2 simulations 
were computed on the domain [—100, 100], and the \i = .1 simulations were computed 
on the domain [-200, 200]. 



http : / /www . math . toronto . edu/simpson/f il.es/media/broadband/model_ 1357 . mp4 



http : / /www . math . toronto . edu/simpson/f iles/media/broadband/mode3_ 1357 .mp4 



http : / /www . math . toronto . edu/simpson/f iles/media/broadband/mode5_1357 . mp4 



http : / /www .math . toronto . edu/simpson/f iles/media/broadband/mode7_ 1357 .mp4 



As one would expect, there is better agreement between the approximate small amplitude 
soliton and the time dependent simulation as jj, — > 0. However, for all values of /i presented, 
there is a persistence of the localization, even if there is distortion to some of the fine structure 
in the higher harmonics. All of this suggests the gap solutions are robust. 

Many other experiments are possible; simulating with more modes, simulating with larger 
values of fi, and seeding the initial conditions of a smaller system into a larger system. In the 



previous work 19 , the exact gap soliton (3.9) was used as an initial condition for successively 



larger truncations of the extended coupled mode system (1.10). 



6 Open Problems 



We conclude this work with a discussion of open problems concerning the existence of non- 
trivial localized solutions of xNLCME and xNLS, arising for the case of a refractive index 
composed of an infinite array of Dirac delta functions. Some of the challenges include: 



Prove the existence of a non-trivial critical point to h, (4.10), the single parameter 
Rayleigh-Ritz approximation, 



Prove the existence of a non-trivial solution to Hq^ uss , (4.2), the Gaussian Rayleigh- 
Ritz approximation, 



Prove the existence of a non-trivial solution to the coupled NLS equations (3.17) 



Prove the existence of a non-trivial solution to the coupled mode equations (3.3). 



By "non-trivial", we mean a solution in which all modes are active and non-zero. It would 
also be of interest to obtain proofs of existence for arbitrarily large finite truncations of these 
problems. Intimately connected with the last two challenges is the question of appropriate 
function spaces. As discussed in Section 4.3, our variational approximations live in the 
function space X s for 1 < s < j ~ 1.26 for which our Theorems [T] ancj2] are stated. The 
upper value on s that ensures that the interval is nonempty is only approximated numerically 
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Figure 7: Surfaces generated from simulations of the coupled mode system (1.10) truncated 



to four modes with initial data (5.1). The \i 
[—50, 50], and the \i 



A simulations were computed on the domain 
.1 simulations were computed on the domain [—200,200]. In both 



cases, there were 2048 grid points. 
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from the "rough" variational approximation. Of course, it is also possible that such solutions 
may not exist. A counterexample would also be of interest. 

Modeling the nonlinear Maxwell equation with refractive index given by a periodic se- 
quence of Dirac delta-functions is a challenging problem both analytically and numerically. 
Results of our work give a starting point to further exploration of this system, and the evo- 
lution of its localized excitations. The question of localized solutions for xNLCME for less 
restrictive, and more physical, refractive indices is also of great interest. 



References 

A. Aceves and S. Wabnitz. Self- induced transparency solitons in nonlinear refractive 
periodic media. Physics Letters A, 141(1-2), 1989. 

D. Agueev and D. Pelinovsky. Modeling of wave resonances in low-contrast photonic 
crystals. SI AM Journal on Applied Mathematics, 65(4): 1101-1129, 2005. 

R. Boyd. Nonlinear Optics. Academic Press, 2008. 

D. Christodoulides and R. Joseph. Slow bragg solitons in nonlinear periodic structures. 
Physical Review Letters, 62(15): 1746-1749, Jan 1989. 

M. Chugunova and D. E. Pelinovsky. Block-diagonalization of the symmetric first-order 
coupled-mode system. SIAM J. Appl. Dyn. Syst, 5(l):66-83, 2006. 

C. de Sterke and J. Sipe. Gap solitons. Progress in Optics, 33:205-205, 1994. 

E. Doedel. AUTO: A program for the automatic bifurcation analysis of autonomous 
systems. Congressus Numerantium, 30:265-284, Jan 1981. 



E. Doedel, B. Oldeman, and et al. AUTO-07P. |http : / / indy . cs . concordia . ca/ auto/ 



B. Eggleton, C. de Sterke, and R. Slusher. Nonlinear propagation in superstructure 
bragg gratings. Optics letters, 21(16):1223-1225, Jan 1996. 

B. Eggleton and R. Slusher. Nonlinear Photonic Crystals. Springer, 2010. 

R. Goodman, R. Slusher, and M. I. Weinstein. Stopping light on a defect. Journal of 
the Optical Society of America B: Optical Physics, 19(7):1635-1652, Jan 2002. 

R. H. Goodman, M. I. Weinstein, and P. Holmes. Nonlinear propagation of light in 
one-dimensional periodic structures. Journal of Nonlinear Science, 11:123-168, 2001. 

D. Pelinovsky and G. Schneider. Justification of the coupled-mode approximation for 
a nonlinear elliptic problem with a periodic potential. Applicable Analysis, 86 (8): 1017- 
1036, 2007. 



31 



[14] D. Pelinovsky and G. Schneider. Moving gap solitons in periodic potentials. Mathemat- 
ical Methods in the Applied Sciences, 31(14): 1739-1760, 2008. 

[15] J. Ranka, R. Windeler, and A. Stentz. Visible continuum generation in air-silica mi- 
crostructure optical fibers with anomalous dispersion at 800 nm. Optics letters, 25(1) :25- 
27, Jan 2000. 

[16] G. Schneider. Nonlinear coupled mode dynamics in hyperbolic and parabolic periodi- 
cally structured spatially extended systems. Asymptotic Analysis, 28 (2): 163-180, 2001. 

[17] G. Schneider and H. Uecker. Existence and stability of modulating pulse solutions in 
Maxwell's equations describing nonlinear optics. Zeitschrift fur Angewandte Mathematik 
und Physik (ZAMP), 54(4):677-712, 2003. 

[18] L. Shampine, I. Gladwell, and S. Thompson. Solving ODEs with MATLAB. Cambridge 
University Press, 2003. 

[19] G. Simpson and M. Weinstein. Coherent structures and carrier shocks in the nonlinear 
maxwell equations. MMS, To Appear. 

[20] R. Tasgal, Y. Band, and B. Malomed. Gap solitons in a medium with third-harmonic 
generation. Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 72(1):1- 
10, Jan 2005. 



32 



